Bioinformática y Estadística 1: Proyecto

# HOMEPAGE

Regresión Bayesiana

Se desea comparar la metodología de regresión lineal simple clásica (frecuentista) con la métodología de regresión Bayeisana. Para esto, considere el siguiente ejemplo.

Ejemplo. Tasa de flujo máximo

La tasa de flujo máximo de una persona es la tasa más rápida a la que la persona puede expulsar aire después de respirar profundamente. El caudal máximo se mide en unidades de litros por minuto y da una indicación de la salud respiratoria de la persona. Para estudiar la relación entre la tasa de flujo máximo y la altura de un individuo, un grupo de investigadores hizo mediciones en una muestra de 17 hombres. Los resultados se muestran en la siguiente tabla. Asuma que la desviación estándar de la variable respuesta es \(\sigma = 115\).

Altura (cm)	Tasa de flujo máximo (l/min)
174	733
183	572
176	500
169	738
183	616
186	787
178	866
175	670
172	550
179	660
171	575
184	577
200	783
195	625
176	470
176	642
190	856

Objetivos

Utilizando el contexto anterior, realice las siguientes actividades.

1. Regresión lineal simple clásica

   1. Ajuste un modelo de regresión lineal simple para predecir la tasa de flujo máximo a partir de la altura.

   2. Muestre la gráfica de la recta de regresión.

   ## Warning in predict.lm(fit, interval = "prediction", level = 0.95): predictions on current data refer to _future_ responses

   

   3. Interprete los parámetros esimtados del modelo usando tanto estimadores puntuales cómo estimadores por intervalo.

   ## 
   ## Call:
   ## lm(formula = Flow ~ Altura, data = D)
   ## 
   ## Residuals:
   ##      Min       1Q   Median       3Q      Max 
   ## -170.096  -99.188    1.904  101.789  216.881 
   ## 
   ## Coefficients:
   ##             Estimate Std. Error t value Pr(>|t|)
   ## (Intercept) -153.926    607.464  -0.253    0.803
   ## Altura         4.511      3.364   1.341    0.200
   ## 
   ## Residual standard error: 115.2 on 15 degrees of freedom
   ## Multiple R-squared:  0.1071, Adjusted R-squared:  0.04757 
   ## F-statistic: 1.799 on 1 and 15 DF,  p-value: 0.1998

El intercepto estimado es \(\approx\) -154 el cual representa el flujo maximo esperado para una persona con altura 0. Este valor por si solo no tiene un sentido real, ya que no existen personas con altura 0. El intercepto actua como un punto de referencia para la recta de regresión.

El coeficiente de la variable explicativa es \(\approx\) 4.5, lo que nos sugiere que por cada incremento de 1cm en la altura de una persona, se espera un incremento de 4.5 l/min en la tasa de flujo máximo. Sin embargo este valor no es significativo, ya que el p-valor es mayor a 0.05.

##                    2.5 %    97.5 %
## (Intercept) -1448.705327 1140.8527
## Altura         -2.657717   11.6807

Dado que incluyen el 0, no son estadisticamente significativos.

Al eliminar el intercepto de la formula se estiman mejor los parametros.

## 
## Call:
## lm(formula = Flow ~ Altura - 1, data = D)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -174.178  -88.720   -2.178   96.142  214.502 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## Altura   3.6601     0.1501   24.39  4.4e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 111.7 on 16 degrees of freedom
## Multiple R-squared:  0.9738, Adjusted R-squared:  0.9722 
## F-statistic:   595 on 1 and 16 DF,  p-value: 4.396e-14

Esto se puede interpretar como que el flujo maximo esperado para una persona con altura 0 es 0, y por cada cm de incremento en la altura se espera un incremento de 3.6 l/min en la tasa de flujo máximo.

##           2.5 %   97.5 %
## Altura 3.342002 3.978201

Estos intervalos no incluyen el 0, por lo que son estadisticamente significativos.

2. Regresión lineal simple Bayesiana

1. Ajuste un modelo de regresión lineal simple Bayesiano para predecir la tasa de flujo máximo a partir de la altura. Para esto, considere una a priori conjunta para el vector \[ \boldsymbol{\beta} \sim N_2(\boldsymbol{\mu}_0, \boldsymbol{\Sigma}_0), \] donde \[ \boldsymbol{\beta}=(\beta_0, \beta_1)^T, \quad \boldsymbol{\mu}_0 = (0,0)^T, \quad \boldsymbol{\Sigma}_0 = \begin{pmatrix} 1000 & 0 \\ 0 & 1000 \end{pmatrix}. \] Luego, la distribución posterior es \(\boldsymbol{\beta}|\boldsymbol{Y}, \boldsymbol{X}, \sigma \sim N_2(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)\), donde \[ \boldsymbol{\Sigma}_1 = \left(\frac{1}{\sigma^2} \boldsymbol{X}^T \boldsymbol{X} + \boldsymbol{\Sigma}_0^{-1}\right)^{-1},\quad \boldsymbol{\mu}_1 = \left( \frac{1}{\sigma^2} \boldsymbol{X}^T \boldsymbol{Y} + \boldsymbol{\Sigma_0}^{-1} \boldsymbol{\mu} \right) \boldsymbol{\Sigma}_1 . \] Recuerde que \(\boldsymbol{X}\) es la matriz de diseño donde los valores de la primera columna son unos y los valores de la segunda columna están dados por la variable explicativa. Más aún, \(\boldsymbol{Y}\) es el vector de observacions de la variable respuesta.

## 
## Media posterior (mu_n): -36.63122 3.862714

## 
## Intervalo de credibilidad del 95% para B0: -92.88197 16.96037

## 
## Intervalo de credibilidad del 95% para B1: 3.573873 4.173469

2. Muestre la gráfica de la recta de regresión que se obtiene con el MAP y compare con el modelo clásico.

3. Muestre e interpreta las distribuciones a posteriori marginales de los parámetros \(\beta_0\) y \(\beta_1\).

La distribucion a posteriori de \(\beta_0\) apunta a que la media es cercana 0 lo que podria concordar con la eliminacion del intercepto en el primer modelo, mientras que la distribucion a posteriori de \(\beta_1\) apunta a que la media es cercana a 3.6 lo que concuerda con el modelo que incluye el intercepto.

4. Muestre la distribución posterior conjunta de \(\beta_0\) y \(\beta_1\). ¿Qué tanta correlación hay entre ambos parámetros?

## La correlación entre beta_0 y beta_1 es: -0.9953782

Como se puede ver en la grafica los valores estan fuertemente relacionados, lo que sugiere que si uno de los parametros es grande, el otro tambien lo sera. El valor que devuelve el comando cor() es \(\approx\) -1 lo que sugiere una fuerte correlacion negativa entre los parametros.

Grafica 3D

3. Simulación Monte Carlo

## Tasa de aceptación para b0: 0.38014 
## Tasa de aceptación para b1: 0.34676

1. Cree un algoritmo Metropolis para simular de la distribución posterior del vector \((\beta_0, \beta_1)\). Explique brevemente la forma de proponer nuevos pasos de la cadena de Markov, además del burn-in y el lag utilizados para generar una muestra pseudo-independiente de tamaño 1000.

Diagrama de flujo

El algoritmo Metropolis-Hastings utiliza una propuesta para sugerir nuevos valores de \(\beta_0\) y \(beta_1\). En este caso, se actualizan los parametros de manera alternada, es decir, en cada iteracion se actualiza un parametro a la vez. Para \(\beta_0\) se propone un nuevo valor a partir de una distribucion uniforme centrada en el valor actual cuyo rango se encuentra en \([-115, 115]\), mientras que para \(\beta_1\) se propone un nuevo valor a partir de una distribucion uniforme centrada en el valor actual cuyo rango es \([-.7,.7]\). El burn-in(1000, ya que es suficiente para ver estabilidad en la cadena) es el numero de iteraciones que se descartan al inicio de la cadena, mientras que el lag(\(\beta_0lag = 1000, ~\beta_1lag = 900\)) es el numero de iteraciones que se descartan entre cada muestra.

2. Compare la distribución muestral de los parámetros contra la función de densidad posterior obtenida en el punto 2.

La discrepancia se debe a que el algoritmo de montecarlo esta aproximado por simulacion, con màs simulaciones deberian de aproximarse mejor.

3. Utilizando las muestras de \((\beta_0, \beta_1)\), cree bandas de credibilidad para la recta de regresión y bandas de predicción. Compare estas bandas con las que se obtienen con el método clásico.

## Warning in predict.lm(fit, interval = "prediction", level = 0.95): predictions on current data refer to _future_ responses

Las dos gráficas muestran bandas de credibilidad muy similares. Sin embargo, las bandas de predicción son diferentes. En la gráfica de regresión lineal clásica, las bandas están más alejadas de la recta, lo que podría indicar que el modelo tiene menor precisión. Por otro lado, en la gráfica generada con los valores obtenidos del algoritmo MCMC, las bandas de predicción son más estrechas, lo que sugiere una mayor precisión en las predicciones.

4. Preguntas específicas.

2. ¿Cuál sería su tasa de flujo máxima promedio para personas con una altura de 1.6 metros? Utilice un intervalo de credibilidad para la media del 95% para responder.

## Intervalo de credibilidad del 95% para la tasa de flujo a 1.6 metros: 451.3561 a 693.7358

(c)Suponga que una persona mide exactamente 2.1 metros. ¿Cuál sería su tasa de flujo máxima de acuerdo al modelo? Utilice un intervalo de credibilidad predictivo del 95% para responder.

## Intervalo de credibilidad predictivo del 95% para la tasa de flujo a 2.1 metros: 734.2064 a 1068.48

4. A partir de los resultados anteriores, ¿le parece que el modelo de regresión es bueno para predecir la tasa de flujo máximo? Explique.

El modelo parece ser razonablemente bueno, ya que:

Los intervalos de credibilidad predictivos no son excesivamente amplios, lo que indica que el modelo está haciendo predicciones bastante precisas para las nuevas observaciones dentro del rango de las alturas que tiene. Ademas que no esta haciendo predicciones demasiado alejadas de lo que se podria esperar con los datos observados.

Aunque el modelo parece funcionar bien, el algoritmo de MCMC tiene puntos para mejorar como el numero de iteraciones para afinarlo mejor, o incrementar el numero de muestras.